

float workRad;
float workTheta;
float workWinding;

struct component_params {
	int type;
	float r0, z0, armStrength, kscale; 
	vec3 RGBmodifier;
};

struct galaxy_type {
	vec3 center;
	vec3 orientation;
	float arm0, arm1;
	int no_components;
	component_params components[10];
};

galaxy_type galaxy;



float getRadius(vec3 rpos, out vec3 P, out float height)  {
  vec3 p = rpos;
  vec3 o = galaxy.orientation;
  height = dot(p,o);
  P = p - o*height ;
  if (abs(P.x*P.x)+ abs(P.y*P.y)+abs(P.z*P.z)<0.00001)
    return 0.00001;
  return length(P);
}





//	bool is = CVector::intersectSphere(RP.camera.camera-g->params.position, RP.direction, g->params.axis, isp1, isp2);

float getHeightModulation(in float height, in float z0) {
    float h = abs(height/z0);
    float val = 0.0;
    if (h>2.0) {
        height = 0.0;
        return 0.0;
    }

    val = 1.0/cosh(h);
    h = val*val;
    return h;
}

float getRadialIntensity(in float rad, in float r0) {
    float r = exp(-rad/r0);
    return clamp(r - 0.01,0.0,1.0);
  }

float findDifference(in float t1, in float t2) {
  float v1 = abs(t1-t2);
  float v2 = abs(t1-t2-6.28);
  float v3 = abs(t1-t2+6.28);
  float v4 = abs(t1-t2-6.28*2.0);
  float v5 = abs(t1-t2+6.28*2.0);

  float v = min(v1,v2);
  v = min(v,v3);
  v = min(v,v4);
  v = min(v,v5);

  return v;

}



float getTheta(in vec3 P) {
  //rotmat.Mul(p, quatRot);
  //return atan2(*rA, *rB) + componentParams.delta;

	return atan(P.x, P.y);
}


float getWinding(in float rad, in float B, in float N) {
  float r = rad + 0.05;
  float t = atan(exp(-0.25/(0.5*(r)))/B)*2.0*N;
  return t;

}



float getArm(in float rad, in vec3 p, in float disp, in float armStrength) {
  float v = abs(findDifference(workWinding,workTheta + disp))/3.14159;
  return pow(1.0-v,armStrength*15.0);
}

float calculateArmValue(in float rad, in vec3 P, in float armStrength, in float arm1, in float arm2) {

  workRad = sqrt(rad);
  workTheta = -getTheta(P);
  workWinding = getWinding(rad, 0.5, 4.0);

  float v1 = getArm(rad, P, arm1, armStrength);
//  if (params->noArms==1)
 //   return v1;
  float v = max(v1, getArm(rad, P, arm2, armStrength));
 // if (params->noArms==2)
    return v;

/*  v = max(v, getArm(rad, P, params->arm3));
  if (params->noArms==3)
    return v;

  v = max(v, getArm(rad, P, params->arm4));

  return v;*/
}


float getPerlinCloudNoise(in vec3 p, in float rad, in float t, in int N, in float ks, in float pers, in int inv, in float freq) {
  //CVector r = twirl(p, t);
  vec3 r = p;
  float perlinnoise = 0.0;
  float max = 0.0;

  float randShiftX = 12513.2351;
  float randShiftY = 5106.96124;
  float std = 1.0;
  vec3 add = vec3(randShiftX, randShiftY, randShiftX*1.521);

  for (int ii=1;ii<N;ii++) {
    float k= (float(ii)*ks);
    float pp = abs(snoise( vec3((r.x + add.x)*k, (r.y+ add.y)*k, (r.z + add.z)*k)));
    add = add + vec3(k*0.111,k*k,-0.712*k);
    if (inv==1) {
      perlinnoise+=(1.0/(pp/k+0.015));
      max+=1.0/(std/k + 0.015);
    }
    else {
      perlinnoise+=pp/k;
      max+=std/k;
    }
  }

  return perlinnoise/max;
  }




float calculateIntensity(in float I, vec3 pos, in component_params c, in int color) {
	float z;
	vec3 P;
	float currentRadius = getRadius(pos, P, z);
	z = getHeightModulation(z, c.z0);
	int N = 20;
    if (z>0.001) 
    {
        float intensity = getRadialIntensity(currentRadius, c.r0);
        if (intensity>0.01) {
        float scale = 1.0;
        float strength = 0.2;
        float armVal = calculateArmValue(currentRadius, P, c.armStrength, galaxy.arm0, galaxy.arm1 );
        float color = c.RGBmodifier[color];
        float val = (strength)*scale*armVal*z*intensity*color;
        float noise = 1.0;
        if (c.type==1) { // disk
        	noise = getPerlinCloudNoise(pos, currentRadius, 0.0, N, c.kscale, 0.0, 1, 0.0);
            return I+val*noise;

        }
        if (c.type==3) { // disk
        	noise = getPerlinCloudNoise(pos, currentRadius, 0.0, N, c.kscale, 0.0, 0, 0.0);
        	noise = pow(noise, 2.0);
        	float s = 250.0;
            return I*exp(-val*noise*s);

        }
        }
    }
    return I;
}

float Integrate(vec3 ray, vec3 from, vec3 to, in int color) {

	vec3 pos = from;
	float I = 0.0;
	float dx = 0.002;
	float l = length(from - to);
	float N = l/dx;
	vec3 d = -normalize(from-to)*dx;
	
	for (int i=0;i<int(N);i++) {
		pos += d;
		for (int j=0;j<galaxy.no_components;j++) {
			I=calculateIntensity(I, pos, galaxy.components[j], color);
		}
	}
	return I;
}



float calculate(vec3 origin, vec3 ray, int color) {

	vec3 isp1, isp2;
	float is = intersectSphere(camera - galaxy.center, ray, vec3(1.0,1.0,1.0), isp1, isp2);
	if (is==1.0) {
		return Integrate(ray, isp1, isp2, color);
			
	}
	return 0.0;
}


void setupGalaxy() {
	galaxy.center = vec3(0,0,0);
	galaxy.orientation = normalize(vec3(0.0,0.0,1.0));
	galaxy.arm0 = 0.0;
	galaxy.arm1 = 3.1415;
	galaxy.no_components = 2;
	
	galaxy.components[0].type = 1;
	galaxy.components[0].r0 = 0.20;
	galaxy.components[0].z0 = 0.02;
	galaxy.components[0].kscale = 1.0;
	galaxy.components[0].armStrength = 0.5;
	galaxy.components[0].RGBmodifier = vec3(1.0, 0.5, 0.2);

	galaxy.components[1].type = 3;
	galaxy.components[1].r0 = 0.25;
	galaxy.components[1].z0 = 0.03;
	galaxy.components[1].kscale = 1.5;
	galaxy.components[1].armStrength = 1.0;
	galaxy.components[1].RGBmodifier = vec3(0.2, 0.3, 1.0);

	
}


void main(void)
{
	setupViewMatrix();
	setupGalaxy();
	vec3 ray = coord2ray(gl_TexCoord[0].x, gl_TexCoord[0].y ,1.0, 1.0, 60.0);
	
	float R = calculate(camera, ray, 0);//Integrate(ray);
	float G = calculate(camera, ray, 1);//Integrate(ray);
	float B = calculate(camera, ray, 2);//Integrate(ray);
	
	vec4 c = vec4(R,G,B,1.0);
//	vec4 c = vec4(camera.x, camera.y, camera.z,1.0);
	gl_FragColor = c;
	
}
